home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / join.pro < prev    next >
Text File  |  1997-07-08  |  10KB  |  408 lines
  1. ; $Id: join.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7.  
  8.  Function EuclidRule, Case1, RCases
  9. ; Dist returns the vector of Euclidean Distances from Case 
  10. ; to the other cases in the array RCase
  11.  
  12.  SC=Size(RCases)
  13.  C=SC(1)
  14.  if(SC(0) GT 1) THEN  R=SC(2) ELSE R=1
  15.  M1= Case1#(Fltarr(R)+1)
  16.  M2=Fltarr(C)+1
  17.  M1=M1-RCases
  18.  if(R EQ 1) THEN  Return,SQRT(Total(M1*M1)) $
  19.   else Return, transpose (SQRT(M2#(M1*M1)))
  20.  END
  21.  
  22. Function Distance1,Case1,RCases,N
  23. ; Determine appropriate rule for computing distances between
  24. ; cases and return the vector of distances between
  25. ; RCases and Case1.
  26.  
  27. SC=Size(RCases)
  28. S=SC(1)
  29. if(SC(0) GE 2) THEN R=SC(2) else R=1
  30.  
  31. Case N of
  32.     "%":BEGIN
  33.        X=Case1#replicate(1.,R) - RCases
  34.        A=where(X NE 0,count)
  35.        if(count NE 0) THEN X(A)=1
  36.        if(SC(0) GT 1) THEN return,       $
  37.          transpose(replicate(1.,S)#X*1/S)    $
  38.        ELSE return,Total(X)*1/S
  39.        END
  40.     
  41.     "COR":BEGIN
  42.           V= [Correlate(Case1,RCases(*,0))]
  43.           for i=1,R-1 DO V=[V,Correlate(Case1,RCases(*,i))]
  44.           return,1-V
  45.           END
  46.     ELSE: return,EuclidRule(Case1,RCases)
  47. ENDCASE
  48. END
  49.  
  50.  
  51.    
  52.        
  53.  
  54.  
  55. Function Normal1, Data,R,C
  56. ;Normal returns Data normalized by columns
  57. Y= Data-Data#Replicate(1./R,R) #Replicate(1.,R)
  58. std =sqrt(Y^2 # Replicate(1./(R-1),R))
  59. D1=Fltarr(c,c)
  60. for i=0,C-1 do     $
  61.   if std(i) NE 0 then    $
  62.      D1(i,i)=1./std(i) else d1(i,i)=0
  63. return, D1#Y
  64. end
  65.  
  66.  
  67. Function FindRow,I,R
  68. ; D is a symmetric matrix stored linearly. FindRow
  69. ; reconstructs the ith row. R is the total number of rows.
  70.  
  71. Common JBlock,D,INDEX
  72.   case I of 
  73.  
  74.   R-1:Begin               ;RETURN,(R-1)*(R-2)/2:(R+1)*(R-2)/2)
  75.       X=LINDGEN((R+1)*(R-2)/2 +1)
  76.       Return,X((R-1)*(R-2)/2:(R+1)*(R-2)/2)
  77.       END
  78.  
  79.   0: BEGIN
  80.      T=LIndgen(R-1)
  81.      RETURN,(1+T)*T/2
  82.      END   
  83.  
  84.   ELSE: BEGIN
  85.         A1=I*(I-1)/2
  86.         A2=(I+1)*(I+2)/2 -1
  87.         T=LIndGen(R-I-1)
  88.         T=(T+1)*T/2
  89.         RETURN,[A1+LIndgen(I),A2 + I*LIndgen(R-I-1 )+T]
  90.         END
  91.  ENDCASE
  92.  END
  93.  
  94.  
  95.  
  96.  Pro MinDist,R,DMin,IMin,I    
  97.  ; MinDist computes the minimum distance from case I to any
  98.  ; other case  using the symmetric distance matrix D. This
  99.  ; distance and the corresponding row number are returned
  100.  ; in DMin and
  101.  ; IMIN respectively.
  102.  
  103.  Common JBLock,D,INDEX
  104.  
  105.  DMin(I)=min(D(FindRow(INDEX(I),R)))
  106.  if(!C LE INDEX(I)-1) THEN IMIN(I)=INDEX(!C) else  $
  107.     IMIn(I)=INDEX(!C+1)
  108.  RETURN
  109.  END
  110.  
  111.  
  112.  
  113.  Pro SetVal,V,I,R  
  114.  
  115.  Common JBlock,D,INDEX
  116.  
  117.  D(FindRow(INDEX(I),R))=v
  118.  RETURN
  119.  END
  120.  
  121.  
  122.   FUNCTION AmalDist,I1,J1,R,pos,Am,cl,here
  123.   ;AmalDist computes the distances from the cluster formed
  124.   ; from I and J to the other cases and clusters as the
  125.   ; minimum distance between members
  126.  
  127.   Common JBlock,D,INDEX
  128.  
  129.    I = INDEX(I1)
  130.    J = INDEX(J1)
  131.    X1=Fltarr(R)
  132.    Y1=Fltarr(R)
  133.  
  134.    
  135.   here = FindRow(I,R)
  136.   X=D(here)
  137.  
  138.   case I of
  139.  
  140.   R-1:X1=[X,1.e30]
  141.   0: X1=[1.e30,X]
  142.   ELSE:X1=[X(0:I-1),1.e30,X(I:R-2)]
  143.   ENDCASE
  144.  
  145.   Y=D(FindRow(J,R))
  146.  
  147.   case J  of
  148.    R-1:Y1=[Y,1.e30]
  149.    0:  Y1=[1.e30,y]
  150.    ELSE:Y1=[Y(0:J-1),1.e30,Y(J:R-2)]
  151.   ENDCASE        
  152.  
  153.     Ind=Where(Y1 EQ 1.e30)
  154.     Y1(Where(X1 EQ 1.e30))= 1.e30
  155.     X1(Ind)=1.e30
  156.  
  157.   INDEX(cl) = I
  158.   INDEX(I) = cl
  159.   
  160.   case Am of
  161.   "MAX" : V=X1>Y1
  162.   "MEAN": V=(pos(I)*X1+pos(J)*Y1)/(pos(I)+pos(J))
  163.    ELse : V=X1<Y1
  164.  ENDCASE
  165.  
  166.  case I of
  167.   0: Return, V(1:*)
  168.   R-1: Return,V(0:R-2)
  169.   ELSE: Return,[V(0:I-1),V(I+1:*)]
  170.  ENDCASE
  171.  
  172.  END
  173.  
  174.  
  175.  
  176.  
  177.  
  178.  Pro Join, DataArray,pos,Imin,Sim,CaseName = CaseName,       $
  179.            List_Name=LN,Norm=NR,Missing=M, Pair=pr,          $
  180.            Distance=Dx, Amal=AM, NoPrint =NP, Width = LS
  181.  
  182.  Common JBlock, D,Index
  183.  
  184. ;+
  185. ; NAME:
  186. ;    JOIN
  187. ;
  188. ; PURPOSE:
  189. ;    To partition the cases represented by the rows in DataArray into
  190. ;    nested clusters. 
  191. ;
  192. ; CATEGORY:
  193. ;    Statistics.
  194. ;
  195. ; CALLING SEQUENCE: 
  196. ;    JOIN, DataArray, pos, Imin, Sim
  197. ;
  198. ; INPUTS: 
  199. ;    DataArray:    a (C,R) dimensioned array where R is the
  200. ;        number of cases to be partitioned and C
  201. ;        is the number of variables.
  202. ;
  203. ; KEYWORDS
  204. ;     CASENAME:    one dimensional array of R case names.
  205. ;
  206. ;    NORM:    Flag to signal whether to normalize the
  207. ;        variable values in Data. Values normalized only if Norm=1. 
  208. ;
  209. ;      MISSING:    Missing data value. If undefined, assume no missing data
  210. ;
  211. ;     DISTANCE:    Rule for computing case distances.
  212. ;        Options are:
  213. ;        1. "%": Distance = percent of variable values that disagree
  214. ;                   between two cases. 
  215. ;        2. "EUCLID":    This is the default.
  216. ;        3. "COR":    Distance between i and    
  217. ;                               j = 1-rij, where rij= the correlation 
  218. ;                between i and j.
  219. ;        4. "OWN":  DataArray= a symmetric a
  220. ;               distance array supplied by the user.  This array
  221. ;               should be a 1-dim array of the elements in the
  222. ;               distance array that are below the main diagonal.
  223. ;
  224. ;    AMAL:    Rule for computing cluster distances. Options are:
  225. ;        1. "MIN": distance = distance between closest members.
  226. ;        2. "MAX": distance = distance between farthest members.
  227. ;        3.  "MEAN" : distance = average of distances between members.
  228. ;                              
  229. ;      NOPRINT:    A flag, if set, to suppress output to the screen.             
  230. ;
  231. ;    WIDTH:    User supplied tree width in characters.  The default is 60.
  232. ;
  233. ; OUTPUT:
  234. ;    The tree of hierarchical clusters is printed to the screen.  
  235. ;     
  236. ; OPTIONAL OUTPUT PARAMETERS: 
  237. ;    Pos:    One-dimensional array of cases in the in the tree of 
  238. ;        hierarchical clusters. Pos(i) = position of case i in tree.
  239. ;
  240. ;    Imin:    One-dimensional array specifying nesting of clusters. 
  241. ;        Imin(i) = smallest cluster properly containing cluster or
  242. ;        case i.
  243. ;              
  244. ;    Sim:    One-dimensional array of minimum distances.
  245. ;        Sim(i) = the distance between the two clusters joined to
  246. ;        form cluster i. 
  247. ;                         
  248. ; RESTRICTIONS:
  249. ;    None.
  250. ;
  251. ; COMMON BLOCKS:
  252. ;    None.
  253. ;
  254. ; SIDE EFFECTS:
  255. ;    None.
  256. ;
  257. ; PROCEDURE:
  258. ;    Adapted from algorithm in Clustering Algorithms by
  259. ;    Hartigan, Wiley Series in Probablity and Mathematical
  260. ;    Statistics, Chapt.4. Function  kmeans1 implements a
  261. ;    function that given a partition P returns a
  262. ;    partition P' in the same neighborhood with reduced
  263. ;    in group error. Function is called repeatedly until 
  264. ;    it finds a fixed point or local minimum. Kmeans1
  265. ;    recomputes cluster means after each reassignment.
  266. ;    Procedure Kmeans successively finds partitions
  267. ;    with the starting partition for K the final partition
  268. ;    for K-1 with the case farthest from its cluster mean
  269. ;    split off to form a new cluster.
  270. ;-
  271.  
  272.  On_Error,2
  273.  SD=Size(DataArray)
  274.  
  275.  if (N_ELements(LN) EQ 0) THEN unit = -1  $
  276.  else openw,unit,/Get,LN
  277.  
  278.  if (SD(0) NE 2) THEN BEGIN
  279.    printf,unit, "Join- Data Matrix must be 2-dimensional"
  280.    goto,quit0
  281.  ENDIF
  282.  
  283.  Data = DataArray
  284.  
  285.  if (N_ELEMENTS(M) NE 0) THEN  BEGIN
  286.      Data = listwise(Data, M)
  287.      SD = size(Data)
  288.      if N_ELEMENTS(Data) le 1 or SD(0) eq 1 THEN BEGIN
  289.         printf,unit,"join- Halting, too many missing data values."
  290.         goto,quit0
  291.      ENDIF
  292.   Endif
  293.  
  294.  C= SD(1)
  295.  R=SD(2)
  296.  Inf= 1.e+30 
  297.  
  298.  if N_ELEMENTS(NR) NE 0 THEN              $
  299.    Data=Normal1(Data,R,C)               ;normalize data
  300.  
  301.  if(N_Elements(Dx) EQ 0) THEN  Dx="EUCLID"  ;if not specified,
  302.                                             ; euclidean
  303.                                             ;distance computed
  304.  if(N_Elements(Am) EQ 0) THEN  Am ="MIN"   
  305.                      ;default amalgamation
  306.                                             ;  policy
  307.     
  308.  
  309.  if(Dx EQ 'OWN') THEN D=Data else Begin 
  310.                               ;user specifies distances 
  311.                                          ;between cases
  312.  
  313.   V=Distance1(Data(*,1),Data(*,0),Dx)
  314.            ;Distance between first two cases print,"V=",V
  315.  
  316.   D=[V]                    ;initialize symmetric Distance
  317.                                 ;matrix D
  318.                                 ;Note that D is 1-dimensional
  319.  
  320.  for i=2l,R-1 DO BEGIN                ;Construct the rest of D
  321.   V=Distance1(Data(*,i),Data(*,0:I-1),Dx )
  322.   D=[D,V]
  323.  ENDFOR
  324. ENDELSE              
  325.  
  326. Index = Findgen(2*R)                  ;initialize indices
  327.            
  328. Sim=Fltarr(2*R)                ;measure of closeness of cases
  329.                                             ;in clusters
  330. DMin=Fltarr(2*R)                            ;
  331. ;DDMin=DMIN
  332. IMin=Intarr(2*R)                           
  333.                                            
  334.  
  335. for I=0l,R-1 DO BEGin
  336.  MinDist,R,DMin,IMin,I           ;Compute min distance between
  337.  ENDFOR                          ; cases and store distance in
  338.                                 ;DMIN(i) and nearest case to i
  339.                                 ; in IMIN(I)
  340.  
  341. RN=R
  342. Pos=[Intarr(RN) +1,Intarr(RN-1)]    
  343.  
  344.  for i= RN,2*RN-2 DO BEGIN           ; Amalgamate clusters
  345.   DM=Min(DMin(0:R-1),MI)         ; Find closest clusters
  346.   Sim(I)=DM
  347.   MJ=IMin(MI)                              
  348.  
  349.   DMin(MI)= inf & IMin(MI)= I & DMin(MJ)=inf & IMin(MJ)=I
  350.                                      ; and combine
  351.  
  352.  
  353.   V= AmalDist(MI,MJ,RN,pos,Am,i,here)             
  354.  
  355.   D(here) = V
  356.   R=R+1
  357. if (i NE (2*RN-2)) THEN BEGIN  
  358.   MinDist,RN,DMin,IMin,I
  359.   Pos(i)   = Pos(MJ)+Pos(MI)
  360.   Ind=where((DMin(0:R-1) GT  V) AND  $
  361.         (DMIN(0:R-1) NE inf) ,count)
  362.          ;Compare distance to new cluster with min  
  363.                              ; and update    
  364.    if(count NE 0) THEN BEGIN
  365.     DMin(Ind) = V(Ind)        
  366.     IMin(Ind) = I
  367.    ENDIF
  368.  
  369.  
  370.  
  371.    SetVal,inf,MJ,RN 
  372.  
  373.    for J= 0l,R-1 DO BEGIN       ;Recompute min distance where
  374.                                ;IMin = one of the amalgamated
  375.                                           ;clusters    
  376.     if(IMin(j) EQ MJ) OR ( IMin(j) EQ MI) THEN BEGIN
  377.       if(DMin(j) NE inf) THEN  MinDist,RN,DMin,IMin,J          
  378.     ENDIF
  379.   ENDFOR
  380. ENDIF
  381.  ENDFOR                              ; Amalgamation completed
  382.  
  383.  
  384.   Pos=[Intarr(RN) +1,Intarr(RN-1)]    
  385.  for i=0l,2*RN-3 DO Pos(IMIN(I))=Pos(IMIN(I))+Pos(I)
  386.                   ;Pos(I)= # of objects in Ith cluster
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  for i=0l,2*RN-3 DO  BEGIN     ;Pos(i) = position of i in tree
  393.               
  394.   j=2*RN-i-3
  395.   temp=Pos(j)
  396.   Pos(j)= Pos(IMin(j))
  397.   Pos(IMin(j))=Pos(IMin(j))-temp
  398.  ENDFOR
  399.  pos=pos(0:RN-1)-1
  400.  if( N_Elements(NP) EQ 0) THEN maketree,pos,imin,SIM,unit, $
  401.                                 CaseName=CaseName,Width = LS
  402.  
  403.  quit0:
  404.  if(unit NE -1) THEN Free_Lun,unit
  405.  RETURN
  406.  END
  407.  
  408.